function error_value = process_error_2d(mesh, pde, U, error_type)
[gauss_points_ref, gauss_weights_ref] = gauss_lengdre_ref("2d");
switch error_type
    case "L2"
        err_E = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
            J_det = det(Jacobi);
            gauss_weights_loc = abs(J_det) * gauss_weights_ref;
            gauss_points_loc(:,1) = V(1,1) + Jacobi(1,1)*gauss_points_ref(:,1) + Jacobi(1,2)*gauss_points_ref(:,2);
            gauss_points_loc(:,2) = V(1,2) + Jacobi(2,1)*gauss_points_ref(:,1) + Jacobi(2,2)*gauss_points_ref(:,2);
            x_hat = (Jacobi(2,2)*(gauss_points_loc(:,1)-V(1,1)) - Jacobi(1,2)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
            y_hat = (-Jacobi(2,1)*(gauss_points_loc(:,1)-V(1,1)) + Jacobi(1,1)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
            phi = basis_function_ref(x_hat, y_hat, mesh.basis_type);
            uh_E = zeros(9,1);
            for i = 1:mesh.Nlb
                uh_E = uh_E + U(mesh.elem(n,i)) .* phi{i};
            end
            f = (pde.u(gauss_points_loc(:,1),gauss_points_loc(:,2)) - uh_E).^2;
            r = gauss_weights_loc * f;
            err_E = err_E + r;
        end
        error_value = sqrt(err_E);
    case "H1"
        err_Ex = 0;
        err_Ey = 0;
        for n = 1:mesh.Ne
            E = mesh.node(mesh.elem(n,:),:);
            V = E(1:3,:);
            Jacobi = [V(2,:)-V(1,:); V(3,:)-V(1,:)]';
            J_det = det(Jacobi);
            gauss_weights_loc = abs(J_det) * gauss_weights_ref;
            gauss_points_loc(:,1) = V(1,1) + Jacobi(1,1)*gauss_points_ref(:,1) + Jacobi(1,2)*gauss_points_ref(:,2);
            gauss_points_loc(:,2) = V(1,2) + Jacobi(2,1)*gauss_points_ref(:,1) + Jacobi(2,2)*gauss_points_ref(:,2);
            x_hat = (Jacobi(2,2)*(gauss_points_loc(:,1)-V(1,1)) - Jacobi(1,2)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
            y_hat = (-Jacobi(2,1)*(gauss_points_loc(:,1)-V(1,1)) + Jacobi(1,1)*(gauss_points_loc(:,2)-V(1,2)))/J_det;
            dNdx = basis_function_ref(x_hat, y_hat, mesh.basis_type, "dx");
            dNdy = basis_function_ref(x_hat, y_hat, mesh.basis_type, "dy");
            phix = cellfun(@(x,y) (Jacobi(2,2)*x - Jacobi(2,1)*y)/J_det, dNdx, dNdy, "UniformOutput", false);
            phiy = cellfun(@(x,y) (-Jacobi(1,2)*x + Jacobi(1,1)*y)/J_det, dNdx, dNdy, "UniformOutput", false);
            uh_Ex = zeros(9,1);
            uh_Ey = zeros(9,1);
            for i = 1:mesh.Nlb
                uh_Ex = uh_Ex + U(mesh.elem(n,i)) .* phix{i};
                uh_Ey = uh_Ey + U(mesh.elem(n,i)) .* phiy{i};
            end
            f1 = (pde.dudx(gauss_points_loc(:,1),gauss_points_loc(:,2)) - uh_Ex).^2;
            f2 = (pde.dudy(gauss_points_loc(:,1),gauss_points_loc(:,2)) - uh_Ey).^2;
            rx = gauss_weights_loc * f1;
            ry = gauss_weights_loc * f2;
            err_Ex = err_Ex + rx;
            err_Ey = err_Ey + ry;
        end
        error_value = sqrt(err_Ex + err_Ey);
    otherwise
        error("Invalid error type.");
end
end